/*! 
 ***************************************************************************
 * \file singlePhase1D.cpp
 * Single phase in porous media in 1D.
 *
 * \par Description
 * This code solves the 1D single phase equation for pressure:
 * \f$ \frac{\phi \mu c}{k}\frac{\partial p}{\partial t} = \nabla^2 p \f$.
 * This equation is solved in the interval \f$ x \in [0,100 cm] \f$.
 * The data used are: 
 * The initial condition for the temperature is \f$ P = 0 \f$
 * and the boundary conditions:  \f$ P_L=2 atm \f$ and \f$ P_R=1 atm\f$.
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 *  \verbatim 
   % make
   % ./singlePhase1D  \endverbatim
 ***************************************************************************
 *  \author    Luis M. de la Cruz [ Mon Mar 17 08:47:37 GMT 2008 ]
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/CDS.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;
  
// Scalar field in one dimension
typedef TunaArray<double,1>::huge ScalarField1D;

int main() 
{
    double length, dt, tolerance; 
    int num_nodes, frequency;
    double porosity, viscosity, permeability, compresibility;
    double PL, PR, Tmax;

// Reading data from the input file
    std::ifstream input_cin ("inputSingle1D"); 
    input_cin >> length
	      >> num_nodes
	      >> porosity
	      >> permeability
	      >> viscosity
	      >> compresibility
	      >> PL
	      >> PR
	      >> dt
	      >> Tmax
	      >> frequency;
    input_cin.close();

    Tmax += dt;

    double Gamma = permeability / (porosity * viscosity * compresibility);

    std::cout << "\n Gamma = " << Gamma;

    StructuredMesh<Uniform<double, 1> > mesh(length, num_nodes);
    float dx = mesh.getDelta(X);
    mesh.print();
    ScalarField1D p ( mesh.getExtentVolumes() ); 
    ScalarField1D pn( mesh.getExtentNodes()   );


    pn = 1.0; // initial condition on nodes
    p = 1.0; // initial condition on volumes
 
    p(0) = pn(0) = 2.0; // boundary condition

    InOut::writeToFile(pn, 0, "./DataSingle1D/pres.", dx);

    
#ifdef WITH_GNUPLOT
    GNUplot plotter;
    plotter("set grid");
    plotter("p \"./DataSingle1D/pres.00000\" title \"Initial cond.\"w lp");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    int pause; std::cin >> pause;
#endif

    SparseMatrix< Diagonal< double, 1> > A(num_nodes); 
  //    DiagonalMatrix< double, 1> A(num_nodes); 
    ScalarField1D              b(num_nodes);

    ScalarEquation< CDS<double, 1> > single_phase(p, A, b, mesh.getDeltas());
    single_phase.setDeltaTime(dt); 
    single_phase.setGamma(Gamma);
    single_phase.setDirichlet(LEFT_WALL); 
    single_phase.setDirichlet(RIGHT_WALL);
    single_phase.print();

    ofstream error_fp("./DataSingle1D/error"), residual_fp("./DataSingle1D/residual");
    int iteration = 0;
    double t = dt;
    double error = 1.0;
    double residual = 1.0;

    while (t <= Tmax) {
      std::cout << "\n ---- Iteration " << ++iteration 
		<< "\t time = " << t
		<< "\t Res = " << residual
		<< "\t Err = " << error;
      single_phase.calcCoefficients();
      Solver::TDMA1D(single_phase);
      residual = single_phase.calcResidual();
      error = single_phase.calcErrorL1();

      error_fp << iteration << "\t" << error << std::endl;
      residual_fp << iteration << "\t" << residual << std::endl;	
 
      single_phase.update();

      t += dt;

      if ( !(iteration % frequency) ) {
	NumUtils::interpolateToNodes(pn, p);
	InOut::writeToFile(pn, 1 , "./DataSingle1D/pres.",dx);
	
#ifdef WITH_GNUPLOT
	std::cout << "\n >---> Press <1> and then <enter> to continue \n";
	std::cin >> pause;	
	plotter("p \"./DataSingle1D/pres.00000\" title \"Initial cond.\" w lp, \"DataSingle1D/pres.00001\" title \"Num. sol.\" w lp");
#endif
      }
    }

    std::cout << "\n\n NOTE: The data were written in the [./DataSingle1D] directory";

    return 0;
}




